import pandas as pd
import numpy as np
from scipy.integrate import trapz
import importlib
import useful
importlib.reload(useful)
import matplotlib
import matplotlib.pyplot as plt
#parameter to use Latex in matplotlib
matplotlib.rcParams.update(matplotlib.rcParamsDefault)
#a plot label should be easy to read
plt.rcParams.update({'font.size': 14})
import matplotlib.image as mpimg
# Load and display the image
img = mpimg.imread("cospectrumModel.png.png")
plt.imshow(img)
plt.axis("off") # Hide axes
plt.show()
import matplotlib.image as mpimg
# Load and display the image
img = mpimg.imread("flowchart.jpg")
# Create a larger figure with higher DPI
plt.figure(figsize=(12, 8), dpi=200) # Increase DPI and figure size
plt.imshow(img)
plt.axis("off") # Hide axes
plt.show()
df = useful.read_day("apr01.dat")
NB: in dataframe df i don't have time (but can use data number) - so I modified the rotation algorithm not to need it
per hour it's 7 tables of u,v,w,T to which apply rotations (it's the seven heights) [an hour is 2.16e5 rows]
# separate in 10 hours
# separate in 7 heights
# on one height:
# - plot u,v,w (just to start)
# - rotate hourly
# - test the rotation
# - plot u,v,w,T, \sigma^2_u, \sigma^2_w, \sigma^2_uw, \sigma^2_wT #kinetcik energy components and fluxes (u,w)
# second part:
# - spectra and cospectra
# - sigma squared from the integral
# - smoothing spectra
# - nice plotting
# - apply model
df_matrix = useful.split_dataframe(df, 10, 7)
df_matrix[0][0]
row_keys = ["h20","h21","h22","h23","h0","h1","h2","h3","h4","h5"] #hour of measurement
col_keys = ["z1", "z2", "z5", "z10", "z15", "z20", "z30"] #livelli FLOSS
df_matrix_dict = useful.split_dataframe_dict(df, 10, 7, row_keys, col_keys)
df_matrix_dict["h5_z30"]
# in (i,j): f"{row_keys[i]}_{col_keys[j]}" to iterate
# or for key, chunk in df_dict.items():
# print(f"\nDataFrame {key}:")
# print(chunk)
let's take a look at the data
useful.plotHour(df_matrix_dict, hour= "h5", height = "z30")
we need to align (x,y,z) and thus (u,v,w) to perpendicular to the floor and along the local hourly pressure gradient and we do so by applying three rotation with the requests:
df_matrix_dict["h5_z30"] = useful.rotateHour(df_matrix_dict, hour= "h5", height = "z30")
df_matrix_dict["h5_z30"]
useful.plotHour(df_matrix_dict, hour= "h5", height = "z30")
useful.plotHour(df_matrix_dict, hour= "h20", height = "z30")
Let's take a look at the horizontal mean u velocity at each height for the hours h5 and h20. The gradient of this vertical profile will be necessary for the cospectrum budget model.
heights = ["z1", "z2", "z5", "z10", "z15", "z20", "z30"]
hour="h5"
umeans = []
for height in heights:
hour_loc = hour + "_" + height
df_matrix_dict[hour_loc] = useful.rotateHour(df_matrix_dict, hour, height )
u = df_matrix_dict[hour_loc].iloc[:,0].values
umeans.append(np.mean(u))
heightsmeters = np.array([1,2,5,10,15,20,30])
# Creazione del grafico
plt.figure(figsize=(4,3))
plt.scatter(umeans,heightsmeters, marker='o', color='coral')
# Miglioramenti estetici
plt.title(r"$\bar{u}(z)$")
plt.ylabel("Height $[m]$")
plt.xlabel(r"$\bar{u} \, [m/s]$")
plt.grid(True, which='both', linestyle='--', linewidth=0.5)
# Mostra il grafico
plt.show()
from scipy.interpolate import CubicSpline
# Create the cubic spline
cs = CubicSpline(heightsmeters,umeans)
# Generate a smooth set of x values for plotting the spline fit
x_fit = np.linspace(heightsmeters.min(), heightsmeters.max(), 200)
y_fit = cs(x_fit)
plt.figure(figsize=(4,3))
plt.scatter(umeans,heightsmeters, marker='o', color='coral')
plt.title(r"$\bar{u}(z)$")
plt.ylabel("Height $[m]$")
plt.xlabel(r"$\bar{u} \, [m/s]$")
plt.grid(True, which='both', linestyle='--', linewidth=0.5)
plt.plot(y_fit, x_fit, label="Cubic Spline Fit", color='black', linewidth=2)
plt.legend()
plt.show()
from scipy.interpolate import UnivariateSpline
# Create a smoothing spline with a chosen smoothing factor (s)
smoothing_factor = 0.5 # Adjust this for more or less smoothing
spline = UnivariateSpline(heightsmeters,umeans, s=smoothing_factor)
y_fit = spline(x_fit) # Smoothed spline
plt.figure(figsize=(4,3))
plt.scatter(umeans,heightsmeters, marker='o', color='coral')
plt.title(r"$\bar{u}(z)$")
plt.ylabel("Height $[m]$")
plt.xlabel(r"$\bar{u} \, [m/s]$")
plt.grid(True, which='both', linestyle='--', linewidth=0.5)
plt.plot(y_fit, x_fit, label="Smoothed Spline Fit", color='black', linewidth=2)
plt.legend()
plt.show()
heights = ["z1", "z2", "z5", "z10", "z15", "z20", "z30"]
hour="h20"
umeans = []
for height in heights:
hour_loc = hour + "_" + height
df_matrix_dict[hour_loc] = useful.rotateHour(df_matrix_dict, hour, height )
u = df_matrix_dict[hour_loc].iloc[:,0].values
umeans.append(np.mean(u))
from scipy.interpolate import UnivariateSpline
# Create a smoothing spline with a chosen smoothing factor (s)
smoothing_factor = 0.5 # Adjust this for more or less smoothing
spline = UnivariateSpline(heightsmeters,umeans, s=smoothing_factor)
y_fit = spline(x_fit) # Smoothed spline
plt.figure(figsize=(4,3))
plt.scatter(umeans,heightsmeters, marker='o', color='coral')
plt.title(r"$\bar{u}(z)$")
plt.ylabel("Height $[m]$")
plt.xlabel(r"$\bar{u} \, [m/s]$")
plt.grid(True, which='both', linestyle='--', linewidth=0.5)
plt.plot(y_fit, x_fit, label="Smoothed Spline Fit", color='black', linewidth=2)
plt.legend()
plt.show()
it is interesting also to take a look at the potential temperature profile $\theta_v = 0.0098*\bar{T}*\text{height}$
heights = ["z1", "z2", "z5", "z10", "z15", "z20", "z30"]
hour="h5"
Tmeans = []
for height in heights:
hour_loc = hour + "_" + height
df_matrix_dict[hour_loc] = useful.rotateHour(df_matrix_dict, hour, height )
T = df_matrix_dict[hour_loc].iloc[:,3].values
Tmeans.append(np.mean(T))
Tmeans = np.array(Tmeans)
heightsmeters = np.array([1,2,5,10,15,20,30])
# Creazione del grafico
plt.figure(figsize=(4,3))
plt.scatter(Tmeans*0.0098*heightsmeters,heightsmeters, marker='o', color='rebeccapurple')
# Miglioramenti estetici
plt.title(r"$\theta_v(z)$, surface is snow")
plt.ylabel("Height $[m]$")
plt.xlabel(r"$\theta_v(z) \, [K\cdot m]$")
#plt.xlim(270,280)
plt.grid(True, which='both', linestyle='--', linewidth=0.5)
# Mostra il grafico
plt.show()
heights = ["z1", "z2", "z5", "z10", "z15", "z20", "z30"]
hour="h20"
Tmeans = []
for height in heights:
hour_loc = hour + "_" + height
df_matrix_dict[hour_loc] = useful.rotateHour(df_matrix_dict, hour, height )
T = df_matrix_dict[hour_loc].iloc[:,3].values
Tmeans.append(np.mean(T))
Tmeans = np.array(Tmeans)
heightsmeters = np.array([1,2,5,10,15,20,30])
# Creazione del grafico
plt.figure(figsize=(4,3))
plt.scatter(Tmeans*0.0098*heightsmeters,heightsmeters, marker='o', color='rebeccapurple')
# Miglioramenti estetici
plt.title(r"$\theta_v(z)$, surface is snow")
plt.ylabel("Height $[m]$")
plt.xlabel(r"$\theta_v(z) \, [K\cdot m]$")
#plt.xlim(270,280)
plt.grid(True, which='both', linestyle='--', linewidth=0.5)
# Mostra il grafico
plt.show()
Before calculating spectra and cospectra, one can evaluate sigmas and covariances, which represents turbulent energies and turbulent fluxes, which are terms of the turbulent kinetic energy equation.
covmat = np.cov(df_matrix_dict["h5_z30"].values,rowvar=False)
print(covmat)
from IPython.display import display, Markdown
a = 1.23
b = 4.56
c = 7.89
d = 0.12
table_latex = f"""
\\
\\begin{{array}}{{|c|c|}}
\\hline
\\textbf{{energia/flussi}} & \\textbf{{valori}} \\\\
\\hline
\\text{{$\sigma^2$$_u$$_u$}} & {covmat[0,0]:.4f} \\\\
\\hline
\\text{{$\sigma^2$$_w$$_w$}} & {covmat[2,2]:.4f} \\\\
\\hline
\\text{{$\sigma^2$$_u$$_w$}} & {covmat[0,2]:.4f} \\\\
\\hline
\\text{{$\sigma^2$$_w$$_T$}} & {covmat[2,3]:.4f} \\\\
\\hline
\\end{{array}}
\
"""
display(Markdown(table_latex))
heights = ["z1", "z2", "z5", "z10", "z15", "z20", "z30"]
hour="h20"
suu = []
sww = []
suw = []
swT = []
for height in heights:
hour_loc = hour + "_" + height
df_matrix_dict[hour_loc] = useful.rotateHour(df_matrix_dict, hour, height )
covmat = np.cov(df_matrix_dict[hour_loc].values,rowvar=False)
suu.append(covmat[0,0])
sww.append(covmat[2,2])
suw.append(covmat[0,2])
swT.append(covmat[2,3])
heightsmeters = np.array([1,2,5,10,15,20,30])
# Set up a 2x2 grid of subplots
fig, axs = plt.subplots(2, 2, figsize=(7,6))
# Plot each quantity on a separate subplot
# First subplot for `suu`
axs[0, 0].plot(suu, heightsmeters, marker='o', color='blue')
axs[0, 0].set_xlabel(r"$\sigma^2_{uu} [m^2/s^2]$")
axs[0, 0].set_ylabel("Height $[m]$")
axs[0, 0].grid(True, which='both', linestyle='--', linewidth=0.5)
# Second subplot for `sww`
axs[0, 1].plot(sww, heightsmeters, marker='o', color='green')
axs[0, 1].set_xlabel(r"$\sigma^2_{ww} [m^2/s^2]$")
axs[0, 1].grid(True, which='both', linestyle='--', linewidth=0.5)
# Third subplot for `suw`
axs[1, 0].plot(suw, heightsmeters, marker='o', color='red')
axs[1, 0].set_xlabel(r"$\overline{u'w'} [m^2/s^2]$")
axs[1, 0].set_ylabel("Height $[m]$")
axs[1, 0].grid(True, which='both', linestyle='--', linewidth=0.5)
# Fourth subplot for `swT`
axs[1, 1].plot(swT, heightsmeters, marker='o', color='purple')
axs[1, 1].set_xlabel(r"$\overline{w'T'} [m^2/s^2]$")
axs[1, 1].grid(True, which='both', linestyle='--', linewidth=0.5)
# Adjust layout and show plot
plt.tight_layout()
plt.suptitle("Covariance components vs Height", y=1.02, fontsize=14)
plt.show()
heights = ["z1", "z2", "z5", "z10", "z15", "z20", "z30"]
hour="h5"
suu = []
sww = []
suw = []
swT = []
for height in heights:
hour_loc = hour + "_" + height
df_matrix_dict[hour_loc] = useful.rotateHour(df_matrix_dict, hour, height )
covmat = np.cov(df_matrix_dict[hour_loc].values,rowvar=False)
suu.append(covmat[0,0])
sww.append(covmat[2,2])
suw.append(covmat[0,2])
swT.append(covmat[2,3])
heightsmeters = np.array([1,2,5,10,15,20,30])
# Set up a 2x2 grid of subplots
fig, axs = plt.subplots(2, 2, figsize=(7,6))
# Plot each quantity on a separate subplot
# First subplot for `suu`
axs[0, 0].plot(suu, heightsmeters, marker='o', color='blue')
axs[0, 0].set_xlabel(r"$\sigma^2_{uu} [m^2/s^2]$")
axs[0, 0].set_ylabel("Height $[m]$")
axs[0, 0].grid(True, which='both', linestyle='--', linewidth=0.5)
# Second subplot for `sww`
axs[0, 1].plot(sww, heightsmeters, marker='o', color='green')
axs[0, 1].set_xlabel(r"$\sigma^2_{ww} [m^2/s^2]$")
axs[0, 1].grid(True, which='both', linestyle='--', linewidth=0.5)
# Third subplot for `suw`
axs[1, 0].plot(suw, heightsmeters, marker='o', color='red')
axs[1, 0].set_xlabel(r"$\overline{u'w'} [m^2/s^2]$")
axs[1, 0].set_ylabel("Height $[m]$")
axs[1, 0].grid(True, which='both', linestyle='--', linewidth=0.5)
# Fourth subplot for `swT`
axs[1, 1].plot(swT, heightsmeters, marker='o', color='purple')
axs[1, 1].set_xlabel(r"$\overline{w'T'} [m^2/s^2]$")
axs[1, 1].grid(True, which='both', linestyle='--', linewidth=0.5)
# Adjust layout and show plot
plt.tight_layout()
plt.suptitle("Covariance components vs Height", y=1.02, fontsize=14)
plt.show()
I need an algorithm that
Then go on spectra/cospectra
Questions:
useful.plotPowerSpectra(df_matrix_dict, hour= "h5", height = "z30", normalize = True,)
useful.plotPowerSpectraPremultiplied_Scatter(df_matrix_dict, hour= "h5", height = "z30", normalize = True)
useful.plotPowerSpectraPremultiplied(df_matrix_dict, hour= "h5", height = "z30", normalize = True)
On high frequencies one can see white noise! it is equidistributed along all frequencies, therefore it is a line with $m=1$ in the premultiplied spectrum
the premultiplied spectrum is an energy. its "peak":
kolmogorov idea: one puts itself in the "reference frame of turbulence"
since it is the negative part we are interested on (usually one normalizes by the integral/covariance, which is negative)
useful.plotAllCospectra(df_matrix_dict, hour= "h5", height = "z30")
useful.plotAllPremultipliedCospectra(df_matrix_dict, hour= "h5", height = "z30")
useful.plotPowerSpectra(df_matrix_dict, hour= "h5", height = "z30", normalize = True, movingAver=True, window_size=50, high_freq = 2*10**-2)
useful.plotPowerSpectraPremultiplied(df_matrix_dict, hour= "h5", height = "z30", normalize = True, movingAver=True, window_size=50, high_freq = 2* 10**-2)
useful.plotAllCospectra(df_matrix_dict, hour= "h5", height = "z30",movingAver=True, window_size=50, high_freq = 2*10**-2)
useful.plotAllPremultipliedCospectra(df_matrix_dict, hour= "h5", height = "z30",movingAver=True, window_size=50, high_freq = 2*10**-2)
uw is energy extraction, we expect it to be a negative quantity -> since I want to focus on the turbulence part, I will divide the cospectrum by its integral, which is negative for uw
To apply the model we need:
uw, uT, ww -> the whole spectra
(uu) -> epsilon from a fit with the hallmark (-5/3) power law exponent of turbulence in the velocity fields.
In order to normalize spectra one can use covariance or the proper integral -> I did the latter
useful.plotSomeCospectraNegAndPos(df_matrix_dict, hour= "h5", height = "z30",movingAver=True, window_size=50, high_freq = 2*10**-2)
useful.plotSomePremCospectraNegAndPos(df_matrix_dict, hour= "h5", height = "z30",movingAver=True, window_size=50, high_freq = 2*10**-2)
from the f*E_uu plot I can see that a -large- window may be chosen between the frequencies $10^{-1}$ e $10^1$. How many data points are there?
sampling_rate = 60 # 60 Hz
T = 3600 # 1hour duration
N = int(T * sampling_rate) # Number of samples
frequencies = np.fft.fftfreq(N, 1/sampling_rate) #freq 0 is the mean!
np.sum((frequencies>=0.1) & (frequencies<=10))
#create a mask
target_freqs = (frequencies>= 0.1) & (frequencies<= 10)
frequencies[target_freqs].size
#calc E_uu
freqs, powerSpectrumSmooth = useful.CalcPowerSpectrum(df_matrix_dict, hour= "h5", height = "z30", normalize = True,
quantity='u', movingAve=True, windowsize=50, highfreq = 2*10**-2)
fig,ax = plt.subplots(1,1,figsize=(6,3))
ax.loglog(freqs,powerSpectrumSmooth, color="rebeccapurple", marker='o', linestyle='none', markersize=2, )
slope = -5/3
y_intercept = 1 # you can change this if you want a different intercept
x_line = np.linspace(0.1, 30, 100) # line x values
y_line = x_line**slope* y_intercept
ax.plot(x_line, y_line, linestyle='--', color='red', label=r'$y \sim x^{-5/3}$ ')
ax.grid(True)
ax.legend()
ax.set_ylim(1e-7, 1e4)
ax.set_xlabel("frequencies [Hz]")
ax.set_ylabel(r"$E_{uu} [m^2/s^2]$")
#now we fit on moving logarithmically spaced windows
from scipy.stats import linregress
# Reset lists if re-running code
slopes, intercepts, r_values = [], [], []
window_size= 500
step = 1
target_freqs = (freqs>= 0.1) & (freqs<= 10)
freqs_interval=freqs[target_freqs]
spectrum_interval=powerSpectrumSmooth[target_freqs]
for i in range(0, freqs_interval.size - window_size + 1, step):
window_spec = spectrum_interval[i:i + window_size]
window_freqs = freqs_interval[i:i + window_size]
# Perform linear regression on each window
result = linregress(np.log10(window_freqs), np.log10(window_spec))
# Store slope, intercept, and R-squared
slopes.append(result.slope)
intercepts.append(result.intercept)
r_values.append(result.rvalue ** 2) # R-squared
intercepts=np.array(intercepts)
slopes = np.array(slopes)
# Display results
# Create a figure with two subplots for the histograms
fig, axs = plt.subplots(1, 2, figsize=(10, 4))
# Plot histogram for slopes
axs[0].hist(slopes, bins=50, color='blue', edgecolor='black')
axs[0].set_title('Histogram of Slopes')
axs[0].set_xlabel('Slope')
axs[0].set_ylabel('Frequency')
# Plot histogram for intercepts
axs[1].hist(intercepts, bins=50, color='green', edgecolor='black')
axs[1].set_title('Histogram of Intercepts')
axs[1].set_xlabel('Intercept')
axs[1].set_ylabel('Frequency')
# Display the histograms
plt.tight_layout()
plt.show()
I choose the fit with "the best -5/3" (square the values then argmin),
then i build up the line to plot it,
and extract the epsilon
# Display results
# Create a figure with two subplots for the histograms
fig, axs = plt.subplots(1, 2, figsize=(10, 4))
Pslopes = (slopes + 5/3)**2
# Plot histogram for slopes
axs[0].hist(Pslopes, bins=50, color='blue', edgecolor='black')
axs[0].set_title('Histogram of Slopes')
axs[0].set_xlabel('Slope')
axs[0].set_ylabel('Frequency')
# Plot histogram for intercepts
axs[1].hist(intercepts, bins=50, color='green', edgecolor='black')
axs[1].set_title('Histogram of Intercepts')
axs[1].set_xlabel('Intercept')
axs[1].set_ylabel('Frequency')
# Display the histograms
plt.tight_layout()
plt.show()
print(slopes[np.argmin(Pslopes)], intercepts[np.argmin(Pslopes)])
#calc E_uu
freqs, powerSpectrumSmooth = useful.CalcPowerSpectrum(df_matrix_dict, hour= "h5", height = "z30", normalize = True,
quantity='u', movingAve=True, windowsize=100, highfreq = 10**-2)
fig,ax = plt.subplots(1,1,figsize=(6,3))
ax.loglog(freqs,powerSpectrumSmooth, color="rebeccapurple", marker='o', linestyle='none', markersize=2, )
slope = slopes[np.argmin(Pslopes)]
y_intercept = 10**intercepts[np.argmin(Pslopes)] # you can change this if you want a different intercept
x_line = np.linspace(0.1, 30, 100) # line x values
y_line = x_line**slope* y_intercept
ax.plot(x_line, y_line, linestyle='--', color='red', label=r'$y \sim x^{-5/3}$ ')
ax.grid(True)
ax.legend()
ax.set_ylim(1e-7, 1e4)
ax.set_xlabel("frequencies [Hz]")
ax.set_ylabel(r"$E_{uu} [m^2/s^2]$")
the model for epsilon is:
$ E_\varepsilon(f)=A \,\bar{u}^{2/3}\varepsilon^{2/3}f^{-5/3} $
with $A=0.15$ for $u$ and $0.2$ for $v,w$
"10**intercepts[np.argmin(Pslopes)]" is $f^{-5/3}$ prefactor
u = df_matrix_dict["h5_z30"].iloc[:,0].values
np.mean(u)
epsilon = ( 10**intercepts[np.argmin(Pslopes)] /0.15 )**(3/2) / np.mean(u)
epsilon
#calc E_uu
freqs, powerSpectrumSmooth = useful.CalcPowerSpectrum(df_matrix_dict, hour= "h5", height = "z30", normalize = True,
quantity='u', movingAve=True, windowsize=100, highfreq = 10**-2)
fig,ax = plt.subplots(1,1,figsize=(6,3))
ax.loglog(freqs,powerSpectrumSmooth, color="rebeccapurple", marker='o', linestyle='none', markersize=2, )
slope = -5/3
x_line = np.linspace(0.1, 30, 100) # line x values
y_line = x_line**slope* 0.15 * np.mean(u)**(2/3) * epsilon**(2/3)
ax.plot(x_line, y_line, linestyle='--', color='orange', label=r'$E_{uu} = A \,\bar{u}^{2/3} \epsilon^{2/3} f^{-5/3}$ ')
ax.grid(True)
ax.legend()
ax.set_ylim(1e-7, 1e4)
ax.set_xlabel("frequencies [Hz]")
ax.set_ylabel(r"$E_{uu} [m^2/s^2]$")
eps = useful.getEpsilon(df_matrix_dict, hour= "h5", height = "z30", plot=True)
hours = ["h20","h21","h22","h23","h0","h1","h2","h3","h4","h5"] #hour of measurement
heights = ["z1", "z2", "z5", "z10", "z15", "z20", "z30"] #livelli FLOSS
hour="h5"
epsilons = []
for height in heights:
hour_loc = hour + "_" + height
df_matrix_dict[hour_loc] = useful.rotateHour(df_matrix_dict, hour, height )
eps = useful.getEpsilon(df_matrix_dict, hour, height , plot=False)
epsilons.append(eps)
hour="h20"
epsilons = []
for height in heights:
hour_loc = hour + "_" + height
df_matrix_dict[hour_loc] = useful.rotateHour(df_matrix_dict, hour, height )
eps = useful.getEpsilon(df_matrix_dict, hour, height , plot=False)
epsilons.append(eps)
# Dati
data = [
[0.06577826542237629, 0.05073024224842888, 0.045734714888747974],
[0.04129473340960273, 0.037393911210352945, 0.03145324377875913],
[0.029444981954559846, 0.028835444383382078, 0.024259755177192876],
[0.018261855843780227, 0.01830264235970057, 0.014910673064156766],
[0.012478869368254451, 0.012840306696191888, 0.01114578570488523],
[0.008523239621536182, 0.00843650217244339, 0.0067336521957292],
[0.0023345902251810675, 0.0019842642782632513, 0.001726669732759961],
]
# Etichette delle righe e intestazioni delle colonne
row_labels = ["z1", "z2", "z5", "z10", "z15", "z20", "z30"]
column_labels = ["u", "v", "w"]
# Genera il codice LaTeX per la tabella
latex_code = "\\begin{array}{lccc}\n"
latex_code += " & " + " & ".join(column_labels) + " \\\\\n"
latex_code += "\\hline\n"
for row_label, row_data in zip(row_labels, data):
row = f"{row_label} & " + " & ".join(f"{x:.6f}" for x in row_data) + " \\\\\n"
latex_code += row
latex_code += "\\end{array}"
print(latex_code)
from IPython.display import display, Markdown
display(Markdown(f"$$\n{latex_code}\n$$"))
from IPython.display import display, Markdown
display(Markdown(f"$$\n{latex_code}\n$$"))
import matplotlib.image as mpimg
# Load and display the image
img = mpimg.imread("cospectrumModel.png.png")
plt.imshow(img)
plt.axis("off") # Hide axes
plt.show()
notes:
# evaluate the model and save some intermediate steps (e.g. the model in its simplified form)
ks, cospModel, modelNOFut, Ewws, Futs = useful.getCospectrumModel_allHeights(df_matrix_dict, hour= "h20")
# Fuw from data
heights = ["z1", "z2", "z5", "z10", "z15", "z20", "z30"]
Fwus = np.empty(7, dtype=object)
for i in range(7):
freqs, cospectrumSmooth = useful.CalcCospectrum(df_matrix_dict, hour="h20", height = heights[i], quantity1='w',
quantity2='u', movingAve=True, windowsize=50, highfreq = 2*10**-2)
Fwus[i] = cospectrumSmooth
in the plots we will filter out negative values for readability
# function to plot the results
def plotModel(altezza):
valid_mask = ~np.isnan(ks[altezza]) & ~np.isinf(ks[altezza]) & \
~np.isnan(Fwus[altezza]) & ~np.isinf(Fwus[altezza]) & \
~np.isnan(cospModel[altezza]) & ~np.isinf(cospModel[altezza]) & \
~np.isnan(Ewws[altezza]) & ~np.isinf(Ewws[altezza]) & \
~np.isnan(Futs[altezza]) & ~np.isinf(Futs[altezza])
# Filter y-values for the left plot (positive values only)
left_mask_Fwus = Fwus[altezza] < 0
left_mask_cospModel = cospModel[altezza] > 0
left_mask_Ewws = Ewws[altezza] > 0
left_mask_Futs = Futs[altezza] > 0
left_mask_modelNOFut = modelNOFut[altezza] > 0
# Set up a figure for just the left plot
fig, ax1 = plt.subplots(figsize=(10, 4))
# Plot the original values on the left subplot with fine lines (only positive y-values)
ax1.loglog(ks[altezza][valid_mask & left_mask_Fwus], -ks[altezza][valid_mask & left_mask_Fwus] * Fwus[altezza][valid_mask & left_mask_Fwus], color="black", linestyle="-", linewidth=0.8, label="data")
ax1.loglog(ks[altezza][valid_mask & left_mask_cospModel], ks[altezza][valid_mask & left_mask_cospModel] * cospModel[altezza][valid_mask & left_mask_cospModel], linestyle="-", linewidth=0.8, label="model")
ax1.loglog(ks[altezza][valid_mask & left_mask_Ewws], ks[altezza][valid_mask & left_mask_Ewws] * Ewws[altezza][valid_mask & left_mask_Ewws], linestyle="-", linewidth=0.8, label="Eww")
ax1.loglog(ks[altezza][valid_mask & left_mask_Futs], ks[altezza][valid_mask & left_mask_Futs] * Futs[altezza][valid_mask & left_mask_Futs], linestyle="-", linewidth=0.8, label="Fut")
ax1.loglog(ks[altezza][valid_mask & left_mask_modelNOFut], ks[altezza][valid_mask & left_mask_modelNOFut] * modelNOFut[altezza][valid_mask & left_mask_modelNOFut], linestyle="-", linewidth=0.8, label="modelNOfut")
ax1.set_xlabel(r"$k_x \, [m^{-1}]$")
ax1.set_ylabel(r"$k_x \cdot F_{uw} \, [m/s^2]$")
ax1.legend(fontsize=8)
ax1.grid(True, which="both", linestyle="--", linewidth=0.5)
# Display the plot
plt.suptitle(f"Comparison h={heights[altezza]}")
plt.tight_layout(rect=[0, 0, 1, 0.95]) # Adjust layout to fit title
plt.show()
plotModel(0)
plotModel(1)
plotModel(2)
plotModel(3)
plotModel(4)
plotModel(5)
plotModel(6)
The behaviour of the model seems to worsen when F_uT is considered. For this reason we will study now $\tau(k)$ with the two models, complete and simplified. Be aware: the form of $\tau(k)$ we hypothesized in the model ($k^{-2/3}$ for each $k$) is a very simple one.
# Load and display the image
img = mpimg.imread("tauofk.png")
plt.imshow(img)
plt.axis("off") # Hide axes
plt.show()
ks, taus , tausSimple = useful.getTau_allHeights(df_matrix_dict, hour="h20")
# Plotting
plt.figure(figsize=(8,4))
for i in range(7):
# Create a mask to filter out negative values in ks[i] and taus[i]
positive_mask = (taus[i] > 0)
# Plot only positive values
plt.loglog(ks[i][positive_mask], taus[i][positive_mask], linestyle="-", linewidth=1.1, label=f"{heights[i]}")
plt.xlabel("$k_x \, [m^{-1}]$")
plt.ylabel(r"$\tau \, [s]$")
plt.title("Complete model, h20")
plt.legend()
plt.grid(True, which="both", linestyle="--", linewidth=0.5)
plt.show()
# Plotting
plt.figure(figsize=(8,4))
for i in range(7):
# Create a mask to filter out negative values in ks[i] and taus[i]
positive_mask = (tausSimple[i] > 0)
# Plot only positive values
plt.loglog(ks[i][positive_mask], tausSimple[i][positive_mask], linestyle="-", linewidth=1.1, label=f"{heights[i]}")
plt.xlabel("$k_x \, [m^{-1}]$")
plt.ylabel(r"$\tau \, [s]$")
plt.title("Simplified model, h20")
plt.legend()
plt.grid(True, which="both", linestyle="--", linewidth=0.5)
plt.show()
analysis were conducted on the hour 5 too, showing even more troubles with the $F_{uT}$ part. at the moment they are not included for brevity.